//__device__ void GaussEliminationGPU(float A[][3], float b[])
//{
//	int i,j,k,ik;
//	float mik, temp;
//	for(k=0; k<3; k++)
//	{
//		mik=-1.0f;
//		for(i=k; i<3; i++)
//		{
//			if(fabs(A[i][k])>mik)
//			{
//				mik=fabs(A[i][k]);
//				ik=i;
//			}
//		}
//		if (mik<1.0e-10f)
//		{
//			b[0]=b[1]=b[2]=0.0f;
//			return;
//		}
//		for(j=k; j<3; j++)
//		{
//			//swap(A[ik][j], A[k][j]);
//			temp = A[ik][j];
//			A[ik][j]=A[k][j];
//			A[k][j]=temp;
//		}
//		//swap(b[ik], b[k]);
//		temp = b[ik];
//		b[ik] = b[k];
//		b[k] = temp;
//		b[k] /= A[k][k];
//		for(i=3-1; i>=k; i--)
//			A[k][i]/=A[k][k];
//		for(i=k+1; i<3; i++)
//		{
//			b[i] -= A[i][k]*b[k];
//			for (j=3-1; j>k; j--)
//				A[i][j] -= A[i][k]*A[k][j];
//		}
//	}
//	for(i=3-1; i>=0; i--)
//	{
//		for (j=i+1; j<3; j++)
//			b[i] -= A[i][j]*b[j];
//	}
//}
//

//__global__ void Silhouette(unsigned char* surIn, size_t pitch)
//{
//	uint pixelX = blockIdx.x*blockDim.x + threadIdx.x;
//	uint pixelY = blockIdx.y*blockDim.y + threadIdx.y;
//	if(pixelX>=cRenderParam.renderWidth || pixelY>=cRenderParam.renderHeight)			return;
//	//uint pixelData=((pixelX<<PIX_BIT)|pixelY); uint pixelInd = pixelY*cRenderParam.renderWidth+pixelX;
//
//	int offset = pixelY*pitch + 16*pixelX;
//	float* pixel = (float*)(surIn + offset);
//	float classId = pixel[3];
//	
//	float norm[3];
//	float xyz[3] = {pixel[0], pixel[1], pixel[2]};
//
//	if(classId == NOPOINT)	
//	{
//		pixel[0] = pixel[1] = pixel[2] = 1.0f;
//		return;
//	}
//	
//
//
//	float errF = 1000.0f;
//	float deltaPoint = 1000.0f;
//	float fValue[3] = {0.0f, 0.0f, 0.0f};
//	float jacobi[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
//	//float norm[3] = {0.0f, 0.0f, 0.0f};
//	//int step[3] = {0, 0, 0};
//	
//	float plane[4] = {0.0f, 0.0f, 0.0f, 0.0f};
//	float xyzBak[3];
//	//float normBak[3];
//
//	bool isHor;
//	if(IsHorSil(classId))
//		isHor = true;
//	else if(IsVerSil(classId)||IsHorVerSil(classId))
//		isHor = false;
//	else return;
//
//	if(!isHor) return;
//
//
//
//	float screenX = 2.0f*(pixelX+0.5f)/cRenderParam.renderWidth-1.0f;
//	float screenY = 1.0f-2.0f*(pixelY+0.5f)/cRenderParam.renderHeight;
//	if(isHor)
//	{		
//		plane[0] = (screenY*cRenderParam.viewProj[0][3]-cRenderParam.viewProj[0][1]);
//		plane[1] = (screenY*cRenderParam.viewProj[1][3]-cRenderParam.viewProj[1][1]);
//		plane[2] = (screenY*cRenderParam.viewProj[2][3]-cRenderParam.viewProj[2][1]);
//		plane[3] = (screenY*cRenderParam.viewProj[3][3]-cRenderParam.viewProj[3][1]);
//	}
//	else 
//	{
//		plane[0] = (screenX*cRenderParam.viewProj[0][3]-cRenderParam.viewProj[0][0]);
//		plane[1] = (screenX*cRenderParam.viewProj[1][3]-cRenderParam.viewProj[1][0]);
//		plane[2] = (screenX*cRenderParam.viewProj[2][3]-cRenderParam.viewProj[2][0]);
//		plane[3] = (screenX*cRenderParam.viewProj[3][3]-cRenderParam.viewProj[3][0]);
//	}
//
//
//
//
//	for(int its=0; its<10; its++)
//	{
//		xyz[0] -= fValue[0];
//		xyz[1] -= fValue[1];
//		xyz[2] -= fValue[2];
//
//		fValue[0] = EvalGenAlgSurf(xyz[0], xyz[1], xyz[2]);
//		EvalGenAlgSurfDf(xyz[0], xyz[1], xyz[2], jacobi[0]);
//
//		norm[0] = jacobi[0][0];
//		norm[1] = jacobi[0][1];
//		norm[2] = jacobi[0][2];
//
//		fValue[1] = EvalGenAlgSurfPolarEye(xyz[0], xyz[1], xyz[2], cRenderParam.eye);
//		jacobi[1][0] = EvalGenAlgSurfPolarEye(xyz[0]+DERIVATIVE_DELTA, xyz[1], xyz[2], cRenderParam.eye);
//		jacobi[1][1] = EvalGenAlgSurfPolarEye(xyz[0], xyz[1]+DERIVATIVE_DELTA, xyz[2], cRenderParam.eye);
//		jacobi[1][2] = EvalGenAlgSurfPolarEye(xyz[0], xyz[1], xyz[2]+DERIVATIVE_DELTA, cRenderParam.eye);
//		jacobi[1][0] = (jacobi[1][0]-fValue[1])*DERIVATIVE_DELTA_INV;
//		jacobi[1][1] = (jacobi[1][1]-fValue[1])*DERIVATIVE_DELTA_INV;
//		jacobi[1][2] = (jacobi[1][2]-fValue[1])*DERIVATIVE_DELTA_INV;
//
//		fValue[2] = plane[0]*xyz[0] + plane[1]*xyz[1] + plane[2]*xyz[2] + plane[3];
//
//		jacobi[2][0] = plane[0];
//		jacobi[2][1] = plane[1];
//		jacobi[2][2] = plane[2];
//
//		errF = fValue[0]*fValue[0]+fValue[1]*fValue[1]+fValue[2]*fValue[2];
//
//		if (errF <= ERROR_FUNC_EPS)					break;
//
//		GaussEliminationGPU(jacobi, fValue);
//
//
//		deltaPoint = fabs(fValue[0])+fabs(fValue[1])+fabs(fValue[2]);
//		if (deltaPoint <= IntervalEPS)			break;	
//	}
//
//	float worldW	   = xyz[0]*cRenderParam.viewProj[0][3] + xyz[1]*cRenderParam.viewProj[1][3] + xyz[2]*cRenderParam.viewProj[2][3] + cRenderParam.viewProj[3][3];
//	fValue[0] = (xyz[0]*cRenderParam.viewProj[0][0] + xyz[1]*cRenderParam.viewProj[1][0] + xyz[2]*cRenderParam.viewProj[2][0] + cRenderParam.viewProj[3][0])/worldW; 
//	fValue[1] = (xyz[0]*cRenderParam.viewProj[0][1] + xyz[1]*cRenderParam.viewProj[1][1] + xyz[2]*cRenderParam.viewProj[2][1] + cRenderParam.viewProj[3][1])/worldW; 
//	fValue[2] = (xyz[0]*cRenderParam.viewProj[0][2] + xyz[1]*cRenderParam.viewProj[1][2] + xyz[2]*cRenderParam.viewProj[2][2] + cRenderParam.viewProj[3][2])/worldW; 
//
//	fValue[0] = ((1.0f+fValue[0])*cRenderParam.renderWidth )/2.0f-0.5f;
//	fValue[1] = ((1.0f-fValue[1])*cRenderParam.renderHeight)/2.0f-0.5f;
//	//float pixelRefine;
//	//if(isHor)		pixelRefine = (fValue[0]-pixelX);
//	//else			pixelRefine = (fValue[1]-pixelY);
//
//	float pixelNew;
//	int pixelOri;
//	int pixelRefine, pixelIter, dir, lim, dir4;
//	if(isHor)
//	{		
//		pixelNew = fValue[0];
//		pixelOri = pixelX;
//		dir = (pixelNew>pixelOri) ? 1 : -1;
//		dir4 = 4*dir;
//		lim = cRenderParam.renderWidth;
//	}
//	else
//	{
//		pixelNew = fValue[1];
//		pixelOri = pixelY;
//		dir = (pixelNew>pixelOri) ? 1 : -1;
//		dir4 = int(pitch/4)*dir;
//		lim = cRenderParam.renderHeight;
//	}
//
//	if(errF<=ERROR_FUNC_EPS||deltaPoint <= IntervalEPS&&abs(pixelNew-pixelOri)<6)
//	{
//		normalize(norm);
//
//		float pad[3], stepN;//, padN[3]
//
//		if(dir>0)	pixelRefine = int(ceil(pixelNew));
//		else		pixelRefine = int(floor(pixelNew));
//		pixelRefine = min(pixelRefine, lim-1);
//		pixelRefine = max(0, pixelRefine);
//
//		pixel-=dir4;
//		//pixelNorm-=dir4;
//		if((pixelX-dir)>=0&&(pixelX-dir)<lim)
//		{
//			stepN = 1.0/(abs(pixelRefine-pixelOri)+1);
//			xyzBak[0] = pixel[0];
//			xyzBak[1] = pixel[1];
//			xyzBak[2] = pixel[2];
//			pad[0] = (xyz[0]-xyzBak[0])*stepN;
//			pad[1] = (xyz[1]-xyzBak[1])*stepN;
//			pad[2] = (xyz[2]-xyzBak[2])*stepN;
//// 			normBak[0] = pixelNorm[0];
//// 			normBak[1] = pixelNorm[1];
//// 			normBak[2] = pixelNorm[2];
//// 			padN[0] = (norm[0]-normBak[0])*stepN;
//// 			padN[1] = (norm[1]-normBak[1])*stepN;
//// 			padN[2] = (norm[2]-normBak[2])*stepN;
//		}
//		else	return;		
//
//		for(pixelIter=pixelOri; (pixelIter!=pixelRefine); pixelIter+=dir, pixel+=dir4/*, pixelNorm+=dir4*/)
//		{
//			xyzBak[0] += pad[0];
//			xyzBak[1] += pad[1];
//			xyzBak[2] += pad[2];
//			pixel[0] = xyzBak[0];
//			pixel[1] = xyzBak[1];
//			pixel[2] = xyzBak[2];
//			pixel[3] = 1.0f;
//			
//// 			normBak[0] += padN[0];
//// 			normBak[1] += padN[1];
//// 			normBak[2] += padN[2];
//// 			pixelNorm[0] = normBak[0];
//// 			pixelNorm[1] = normBak[1];
//// 			pixelNorm[2] = normBak[2];
//		}
//
//		//pixel = (float*)(surIn + offset);
//		pixel[0] = xyz[0];
//		pixel[1] = xyz[1];
//		pixel[2] = xyz[2];
//
//	
//		//OUTPUTVEC(xyz);
//		//pixel[3] = (pixelNew-pixelRefine)*dir;
//
//	
//		//pixelNorm[0] = norm[0];
//		//pixelNorm[1] = norm[1];
//		//pixelNorm[2] = norm[2];
//		//pixelNorm[3] = isHor;
//	}
//
//}
//
//#define KIONE_DONE -1.0
//#define KIONE_SIL_DONE -2.0
//__global__ void Classify(unsigned char* Norm, unsigned char* Pos, size_t pitch, float* debug, float* FillHor, float* FillVer, uint* SilScan, size_t* num)
//{	
//	uint x = blockIdx.x*blockDim.x + threadIdx.x;
//	uint y = blockIdx.y*blockDim.y + threadIdx.y;
//	if(x>=cRenderParam.renderWidth || y>=cRenderParam.renderHeight)			return;
//
//	int offset = y*pitch + 16*x;
//	float* pNorm = (float*)(Norm + offset);
//	// = {pixel[0], pixel[1], pixel[2]};
//	float classIdNorm = pNorm[3];
//
//
//	float* pixel = (float*)(Pos + offset);
//	float classId = pixel[0];
//	
//	float norm[3];
//	float xyz[3] = {pixel[1], pixel[2], pixel[3]};
//
//	//pixel[0] = pixel[1] = pixel[2] = 1.0f;
//
//	if(classId == NOPOINT)	
//	{
//		pNorm[3] = KIONE_DONE;
//		return;
//	}
//
//	if(IsSil(classId))
//	{
//		int ind; uint silInfo=((x<<PIX_BIT)|y)<<1;
//		if(classIdNorm*.8<1)//vertical
//		{
//			ind = atomicAdd(num, 1);
//			SilScan[ind] = silInfo;
//		}
//		if(classIdNorm>0.8)//horizon
//		{
//			ind = atomicAdd(num, 1);
//			SilScan[ind] = silInfo+1;
//		}
//		pNorm[3] = KIONE_DONE;
//		return;		
//	}
//
//	float theta[3];
//	theta[0] = cRenderParam.eye[0]-xyz[0];
//	theta[1] = cRenderParam.eye[1]-xyz[1];
//	theta[2] = cRenderParam.eye[2]-xyz[2];
//	normalize(theta);
//
//
//	if(!NewtonGAS(xyz, norm, theta)) 
//	{
//		return;
//	}
//	
//	float nearSil = theta[0]*norm[0]+theta[1]*norm[1]+theta[2]*norm[2];
//	if((nearSil)<0)
//	{
//		pNorm[3] = KIONE_DONE;
//		return;
//	}
//	normalize(norm);
//
//	MATERIAL material;
//	float rgb[3] = {0.0, 0.0, 0.0};
//	chooseMaterial(0, &material, &cRenderParam.material);
//	shading(rgb, xyz, norm, cRenderParam.eye, cRenderParam.nLight, cRenderParam.light, &material);
//	pixel[0] = rgb[0];
//	pixel[1] = rgb[1];
//	pixel[2] = rgb[2];
//	pixel[3] = KIONE_DONE;
//
//
//	pNorm[0] = xyz[0];
//	pNorm[1] = xyz[1];
//	pNorm[2] = xyz[2];
//	pNorm[3] = KIONE_DONE;
//}
//
//void GetSingScreen(float* singScreen)
//{
//	float x,y,w, x1,y1,z1;//z,
//	int ind=0;
//	for(int i=0; i<GASData.nSingPt; i++)
//	{
//		x1 = GASData.singPt[ind++];
//		y1 = GASData.singPt[ind++];
//		z1 = GASData.singPt[ind++];
//
//		x = x1*renderParam->viewProj[0][0]+y1*renderParam->viewProj[1][0]+z1*renderParam->viewProj[2][0]+renderParam->viewProj[3][0];
//		y = x1*renderParam->viewProj[0][1]+y1*renderParam->viewProj[1][1]+z1*renderParam->viewProj[2][1]+renderParam->viewProj[3][1];
//		//z = x1*renderParam->viewProj[0][2]+y1*renderParam->viewProj[1][2]+z1*renderParam->viewProj[2][2]+renderParam->viewProj[3][2];
//		w = x1*renderParam->viewProj[0][3]+y1*renderParam->viewProj[1][3]+z1*renderParam->viewProj[2][3]+renderParam->viewProj[3][3];
//		
//		x1 = x/w;
//		y1 = y/w;
//
//		singScreen[2*i  ] = (x1+1.0f)*width/2.0f-0.5f;
//		singScreen[2*i+1] = (1.0f-y1)*height/2.0f-0.5f;
//		//printf("%f %f\n", singScreen[2*i  ], singScreen[2*i+1]);
//	}
//}
//
//
//
//
//__global__ void UpdataSil(uint* SilScan, unsigned char* surIn, unsigned char* surOut, size_t pitch, float* FillHor, float* FillVer, uint Number, float* debug)
//{
//	uint tid = blockIdx.x*blockDim.x + threadIdx.x;
//	if(tid>=Number)	return;
//
//	uint silInfo = SilScan[tid];
//	bool isHor = silInfo%2; 
//	int x = silInfo>>(PIX_BIT+1);
//	int y = silInfo<<(32-PIX_BIT-1)>>(32-PIX_BIT);
//	//if(x!=85||y!=220)	return;
//	//if(x!=185||y!=106)	return;
//	
//	float* pixel = (float*)(surOut+y*pitch)+4*x;
//	float xyz[3] = {pixel[1], pixel[2], pixel[3]};
//
//	float screenX = 2.0f*(x+0.5f)/cRenderParam.renderWidth-1.0f;
//	float screenY = 1.0f-2.0f*(y+0.5f)/cRenderParam.renderHeight;
//	float plane[4];
//
//
//	if(isHor)
//	{		
//		plane[0] = (screenX*cRenderParam.viewProj[0][3]-cRenderParam.viewProj[0][0]);
//		plane[1] = (screenX*cRenderParam.viewProj[1][3]-cRenderParam.viewProj[1][0]);
//		plane[2] = (screenX*cRenderParam.viewProj[2][3]-cRenderParam.viewProj[2][0]);
//		plane[3] = (screenX*cRenderParam.viewProj[3][3]-cRenderParam.viewProj[3][0]);
//	}
//	else 
//	{
//		plane[0] = (screenY*cRenderParam.viewProj[0][3]-cRenderParam.viewProj[0][1]);
//		plane[1] = (screenY*cRenderParam.viewProj[1][3]-cRenderParam.viewProj[1][1]);
//		plane[2] = (screenY*cRenderParam.viewProj[2][3]-cRenderParam.viewProj[2][1]);
//		plane[3] = (screenY*cRenderParam.viewProj[3][3]-cRenderParam.viewProj[3][1]);
//	}
//
//	float fValue[3] = {0,0,0}, errF, deltaPoint;
//	float Jacobi[3][3], float norm[3];
//
//	for(int its=0; its<cRenderParam.nIter; its++)
//	{
//		xyz[0] -= fValue[0];
//		xyz[1] -= fValue[1];
//		xyz[2] -= fValue[2];
//
//		fValue[0] = EvalGenAlgSurf(xyz[0], xyz[1], xyz[2]);
//		//Jacobi[0][0] = EvalGenAlgSurf(xyz[0]+DERIVATIVE_DELTA, xyz[1], xyz[2]);
//		//Jacobi[0][1] = EvalGenAlgSurf(xyz[0], xyz[1]+DERIVATIVE_DELTA, xyz[2]);
//		//Jacobi[0][2] = EvalGenAlgSurf(xyz[0], xyz[1], xyz[2]+DERIVATIVE_DELTA);
//		//Jacobi[0][0] = (Jacobi[0][0]-fValue[0])*DERIVATIVE_DELTA_INV;
//		//Jacobi[0][1] = (Jacobi[0][1]-fValue[0])*DERIVATIVE_DELTA_INV;
//		//Jacobi[0][2] = (Jacobi[0][2]-fValue[0])*DERIVATIVE_DELTA_INV;
//		EvalGenAlgSurfDf(xyz[0], xyz[1], xyz[2], Jacobi[0]);
//		norm[0] = Jacobi[0][0];
//		norm[1] = Jacobi[0][1];
//		norm[2] = Jacobi[0][2];
//
//		fValue[1] = EvalGenAlgSurfPolarEye(xyz[0], xyz[1], xyz[2], cRenderParam.eye);
//		Jacobi[1][0] = EvalGenAlgSurfPolarEye(xyz[0]+DERIVATIVE_DELTA, xyz[1], xyz[2], cRenderParam.eye);
//		Jacobi[1][1] = EvalGenAlgSurfPolarEye(xyz[0], xyz[1]+DERIVATIVE_DELTA, xyz[2], cRenderParam.eye);
//		Jacobi[1][2] = EvalGenAlgSurfPolarEye(xyz[0], xyz[1], xyz[2]+DERIVATIVE_DELTA, cRenderParam.eye);
//		Jacobi[1][0] = (Jacobi[1][0]-fValue[1])*DERIVATIVE_DELTA_INV;
//		Jacobi[1][1] = (Jacobi[1][1]-fValue[1])*DERIVATIVE_DELTA_INV;
//		Jacobi[1][2] = (Jacobi[1][2]-fValue[1])*DERIVATIVE_DELTA_INV;
//
//		fValue[2] = plane[0]*xyz[0] + plane[1]*xyz[1] + plane[2]*xyz[2] + plane[3];
//		Jacobi[2][0] = plane[0];
//		Jacobi[2][1] = plane[1];
//		Jacobi[2][2] = plane[2];
//
//		errF = fValue[0]*fValue[0]+fValue[1]*fValue[1]+fValue[2]*fValue[2];
//		if (errF <= cRenderParam.funcErr)					break;
//
//		GaussEliminationGPU(Jacobi, fValue);
//
//
//		deltaPoint = fabs(fValue[0])+fabs(fValue[1])+fabs(fValue[2]);
//		if (deltaPoint <= cRenderParam.deltaEps)			break;	
//	}
//	//if(errF > cRenderParam.funcErr)	return;
//
//
//	normalize(norm);
//	float normSrx = cRenderParam.invViewProj[0][0]*norm[0]+cRenderParam.invViewProj[0][1]*norm[1]+cRenderParam.invViewProj[0][2]*norm[2];
//	float normSry = -(cRenderParam.invViewProj[1][0]*norm[0]+cRenderParam.invViewProj[1][1]*norm[1]+cRenderParam.invViewProj[1][2]*norm[2]);
//
//	//z = xyz[0]*cRenderParam.viewProj[0][2]+xyz[1]*cRenderParam.viewProj[1][2]+xyz[2]*cRenderParam.viewProj[2][2]+cRenderParam.viewProj[3][2];
//	int pOffset, silOffset;
//	float part;
//	float w = xyz[0]*cRenderParam.viewProj[0][3]+xyz[1]*cRenderParam.viewProj[1][3]+xyz[2]*cRenderParam.viewProj[2][3]+cRenderParam.viewProj[3][3];
//	if(isHor)
//	//if(abs(normSrx)<abs(normSry))
//	{
//		//abs(normSry/normSrx)>0.8f
//		//if(abs(normSry/normSrx)<1)return;
//		float ySr = xyz[0]*cRenderParam.viewProj[0][1]+xyz[1]*cRenderParam.viewProj[1][1]+xyz[2]*cRenderParam.viewProj[2][1]+cRenderParam.viewProj[3][1];
//		ySr = ySr/w;
//		if(ySr<-1.0||ySr>1.0)	return;
//		ySr = (1.0f-ySr)*cRenderParam.renderHeight/2.0f; //Attention!!!!!!!!!!!!! -0.5f;
//		
//		part = ySr-truncf(ySr);
//		int ySrI = (int)(truncf(ySr))+1;
//		if(normSry<0)  {part = 1-part;ySrI--;}
//		pOffset = ySrI*pitch+16*x;
//		silOffset = abs(y-ySrI);
//
//		int xoff = x*MAX_FILL_N;
//		int ind = atomicAdd((int*)(FillHor+xoff), 1);
//		if(ind+1<MAX_FILL_N)	FillHor[xoff+1+ind] = ySr;
//	}
//	else
//	{
//		//if(abs(normSry/normSrx)>1)return;
//		float xSr = xyz[0]*cRenderParam.viewProj[0][0]+xyz[1]*cRenderParam.viewProj[1][0]+xyz[2]*cRenderParam.viewProj[2][0]+cRenderParam.viewProj[3][0];
//		xSr = xSr/w;
//		if(xSr<-1.0||xSr>1.0)	return;
//		xSr = (xSr+1.0f)*cRenderParam.renderWidth/2.0f; //Attention!!!!!!!!!!!!! -0.5f;
//		part = xSr;
//
//		part = xSr-truncf(xSr);
//		int xSrI = (int)(truncf(xSr))+1;
//		if(normSrx<0)  {part = 1-part;xSrI--;};
//		pOffset = y*pitch+16*xSrI;
//		silOffset = abs(x-xSrI);
//
//		int yoff = y*MAX_FILL_N;
//		int ind = atomicAdd((int*)(FillVer+yoff), 1);
//		if(ind+1<MAX_FILL_N)	FillVer[yoff+1+ind] = xSr;
//	}
//
//	pixel = (float*)(surOut+pOffset);
//	MATERIAL material;
//	float rgb[3] = {0.0, 0.0, 0.0};
//	chooseMaterial(0, &material, &cRenderParam.material);
//	shading(rgb, xyz, norm, cRenderParam.eye, cRenderParam.nLight, cRenderParam.light, &material);
//	//pixel[0] = pixel[1] = pixel[2] = 1.0f;return;
//	pixel[0] = rgb[0]*part;
//	pixel[1] = rgb[1]*part;
//	pixel[2] = rgb[2]*part;	
//	//pixel[0]=pixel[1]=pixel[2]=1.0f; 
//
//	float* pNorm = (float*)(surIn+pOffset);
//	pNorm[0] = xyz[0];
//	pNorm[1] = xyz[1];
//	pNorm[2] = xyz[2];
//	pNorm[3] = KIONE_SIL_DONE;
//
//
//	pOffset = (abs(normSrx)>abs(normSry))?(normSrx<0?4:-4):(normSry<0?pitch/4:-pitch/4);
//	for(int i=0; i<silOffset; i++)
//	{
//		pixel += pOffset;
//		//if(pixel[3]!=KIONE_DONE)
//		{
//			pixel[0] = pixel[1] = pixel[2] = 1.0f;
//		}
//	}
//
//}
//
//void SilTopo(float* singScreen)
//{
//	int head,tail,branch=0, lastBr=0;
//	float* pData = hFillHor;
//	
//	head = 0;
//	while(*((int*)pData)==0)	{pData+=MAX_FILL_N;head++;}
//	for(int i=head; i<width; i++) //Horizonal silhouette are in x-direction.
//	{
//		branch = *((int*)pData);
//		if(branch != lastBr)     
//		{}
//		lastBr = branch;
//	}
//
//}

__device__ bool NewtonGAS2(float xyz[3], float norm[3], float theta[3])
{
	float fValue[3]={0,0,0},errF,deltaPoint;
	float pol[3] = {xyz[0], xyz[1], xyz[2]};
	
	int its;
	for(its=0; its<cRenderParam.nIter; its++)
	{
		pol[0] -= fValue[0];
		pol[1] -= fValue[1];
		pol[2] -= fValue[2];

		fValue[0] = EvalGenAlgSurfPolarEye(pol[0], pol[1], pol[2], cRenderParam.eye);
		norm[0] = EvalGenAlgSurfPolarEye(pol[0]+DERIVATIVE_DELTA, pol[1], pol[2], cRenderParam.eye);
		norm[1] = EvalGenAlgSurfPolarEye(pol[0], pol[1]+DERIVATIVE_DELTA, pol[2], cRenderParam.eye);
		norm[2] = EvalGenAlgSurfPolarEye(pol[0], pol[1], pol[2]+DERIVATIVE_DELTA, cRenderParam.eye);
		norm[0] = (norm[0]-fValue[0])*DERIVATIVE_DELTA_INV;
		norm[1] = (norm[1]-fValue[0])*DERIVATIVE_DELTA_INV;
		norm[2] = (norm[2]-fValue[0])*DERIVATIVE_DELTA_INV;


		errF = fValue[0]*fValue[0];
		if ( errF <= cRenderParam.funcErr)		break;
		fValue[0] = (fValue[0])/(norm[0]*theta[0]+norm[1]*theta[1]+norm[2]*theta[2]);
		fValue[1] = theta[1]*fValue[0];
		fValue[2] = theta[2]*fValue[0];
		fValue[0] = theta[0]*fValue[0];

		deltaPoint = fabs(fValue[0])+fabs(fValue[1])+fabs(fValue[2]);
		if (deltaPoint <= cRenderParam.deltaEps)			break;	
	}

	if(errF > cRenderParam.funcErr)	return false;
	//else return true;

	xyz[0]+=theta[0];
	xyz[1]+=theta[1];
	xyz[2]+=theta[2];
	//return true;


	float b0=1.0, a0=0.0;
	float dir[3]={pol[0]-xyz[0],pol[1]-xyz[1],pol[2]-xyz[2]};
	float F, G, w1, fw1;
	F = EvalGenAlgSurf(xyz[0], xyz[1], xyz[2]);
	G = EvalGenAlgSurf(xyz[0]+b0*dir[0], xyz[1]+b0*dir[1], xyz[2]+b0*dir[2]);
	//if(F*G>0)	return;
	its=0;
	while (its<cRenderParam.nIter)
	{
		w1 = (G*a0-F*b0)/(G-F);
		fw1 = EvalGenAlgSurf(xyz[0]+w1*dir[0], xyz[1]+w1*dir[1], xyz[2]+w1*dir[2]);
		if(fabs(fw1)<ERROR_FUNC_EPS)	{xyz[0]+=w1*dir[0], xyz[1]+=w1*dir[1], xyz[2]+=w1*dir[2]; break;}
		if(F*fw1<0)
		{
			b0=w1;
			G=fw1;
			F=F/2.0;
		}
		else
		{
			a0=w1;
			F=fw1;
			G=G/2;
		}
		its++;
	}
	EvalGenAlgSurfDf(xyz[0], xyz[1], xyz[2], norm);
	return true;

}

__device__ bool NewtonGAS(float xyz[3], float norm[3], float theta[3])
{
	float fValue[3]={0,0,0},errF,deltaPoint;
	
	for(int its=0; its<cRenderParam.nIter; its++)
	{
		xyz[0] -= fValue[0];
		xyz[1] -= fValue[1];
		xyz[2] -= fValue[2];

		fValue[0] = EvalGenAlgSurf(xyz[0], xyz[1], xyz[2]);
		//norm[0] = EvalGenAlgSurf(xyz[0]+DERIVATIVE_DELTA, xyz[1], xyz[2]);
		//norm[1] = EvalGenAlgSurf(xyz[0], xyz[1]+DERIVATIVE_DELTA, xyz[2]);
		//norm[2] = EvalGenAlgSurf(xyz[0], xyz[1], xyz[2]+DERIVATIVE_DELTA);
		//norm[0] = (norm[0]-fValue[0])*DERIVATIVE_DELTA_INV;
		//norm[1] = (norm[1]-fValue[0])*DERIVATIVE_DELTA_INV;
		//norm[2] = (norm[2]-fValue[0])*DERIVATIVE_DELTA_INV;
		EvalGenAlgSurfDf(xyz[0], xyz[1], xyz[2], norm);


		errF = fValue[0]*fValue[0];
		if ( errF <= cRenderParam.funcErr)		break;
		fValue[0] = (fValue[0])/(norm[0]*theta[0]+norm[1]*theta[1]+norm[2]*theta[2]);
		fValue[1] = theta[1]*fValue[0];
		fValue[2] = theta[2]*fValue[0];
		fValue[0] = theta[0]*fValue[0];

		deltaPoint = fabs(fValue[0])+fabs(fValue[1])+fabs(fValue[2]);
		if (deltaPoint <= cRenderParam.deltaEps)			break;	
	}

	if(errF > cRenderParam.funcErr)	return false;
	else return true;

}

#define PIX_BIT 11


__global__ void FirstNewton(unsigned char* surIn, unsigned char* Normal, unsigned char* Pos, size_t pitch, uint* Valid, uint* Data)
{	
	uint x = blockIdx.x*blockDim.x + threadIdx.x;
	uint y = blockIdx.y*blockDim.y + threadIdx.y;
	if(x>=cRenderParam.renderWidth || y>=cRenderParam.renderHeight)			return;
	uint pixelData=((x<<PIX_BIT)|y); uint pixelInd = y*cRenderParam.renderWidth+x;

	int offset = y*pitch + 16*x;
	float* pixel = (float*)(surIn + offset);
	float classId = pixel[3];
	
	float norm[3];
	float xyz[3] = {pixel[0], pixel[1], pixel[2]};


	if(classId == NOPOINT)	
	{
		Valid[pixelInd] = 1;
		Data[pixelInd] = pixelData;
		pixel[0] = pixel[1] = pixel[2] = 1.0f;
		return;
	}
	//return;
	

	
	float theta[3];
	theta[0] = cRenderParam.eye[0]-xyz[0];
	theta[1] = cRenderParam.eye[1]-xyz[1];
	theta[2] = cRenderParam.eye[2]-xyz[2];
	normalize(theta);



	if(!NewtonGAS2(xyz, norm, theta)) 
	//if(!NewtonGAS(xyz, norm, theta)) 
	{
		Valid[pixelInd] = 1;
		Data[pixelInd] = pixelData;

		return;
	}
	
	normalize(norm);
	float nearSil = theta[0]*norm[0]+theta[1]*norm[1]+theta[2]*norm[2];
	//if((nearSil)*cRenderParam.nIter<1)
	if((nearSil)<0)
	{
		Valid[pixelInd] = 1;
		Data[pixelInd] = pixelData;

		pixel[0] = 1.0;
		pixel[1] = 0.0;
		pixel[2] = 0.0;
		pixel[3] = 1.0;

		return;
	}

	pixel[0] = xyz[0];
	pixel[1] = xyz[1];
	pixel[2] = xyz[2];
	pixel[3] = 1.0;
	
	pixel = (float*)(Normal + y*pitch) + 4*x;
	pixel[0] = norm[0];
	pixel[1] = norm[1];
	pixel[2] = norm[2];

	//float rgb[3] = {0.0, 0.0, 0.0};
	////MATERIAL material;
	////chooseMaterial(0, &material, &cRenderParam.material);
	//shading(rgb, xyz, norm, cRenderParam.eye, cRenderParam.nLight, cRenderParam.light, &cRenderParam.material, 0);
	//pixel[0] = rgb[0]*.8;
	//pixel[1] = rgb[1]*.8;
	//pixel[2] = rgb[2]*.8;
}




void compact(uint* compactNum, uint allNum, size_t* dNum, uint* Occupied, uint* OccupiedValid, uint* Compact, uint saveSize, uint* SaveData)
{
	CUDPPConfiguration config;
	CUDPPResult result = CUDPP_SUCCESS;	

	CUDA_SAFE_CALL( cudaMemset(Compact, 0, allNum) );
	
	config.algorithm = CUDPP_COMPACT;
    config.datatype = CUDPP_UINT;
	config.options = CUDPP_OPTION_FORWARD;
	result = cudppPlan(&plan, config, allNum/sizeof(uint), 1, 0);
	if (result != CUDPP_SUCCESS)
	{
		printf("Error creating plan for compact\n");
		return;
	}
	result = cudppCompact(plan, Compact, dNum, Occupied, OccupiedValid, allNum/sizeof(uint));
	if (result != CUDPP_SUCCESS)
	{
		printf("Error compacting !\n");
		return;
	}	
	result = cudppDestroyPlan(plan);    
    if (result != CUDPP_SUCCESS)
    {
        printf("Error destroying CUDPPPlan!\n");
		return;
    }

	CUDA_SAFE_CALL( cudaMemcpy(compactNum, dNum, sizeof(size_t), cudaMemcpyDeviceToHost) );

#if DEBUG_BUFFERS
	printf("dOccupiedCompact:\n");
	dumpBuffer(Compact, *compactNum);
#endif

	if(saveSize!=0)
	{
		if(saveSize<(*compactNum))
			printf("Malloc more memory for SaveData!\n");

		CUDA_SAFE_CALL( cudaMemset(SaveData, 0, saveSize*sizeof(uint)) );
		CUDA_SAFE_CALL( cudaMemcpy(SaveData, Compact, (*compactNum)*sizeof(uint), cudaMemcpyDeviceToDevice) );
		CUDA_SAFE_CALL( cudaMemset(Compact, 0, allNum) );
	}
}


__global__ void SecondKI(uint allNum, uint* Data, unsigned char* Normal, unsigned char* surface, size_t pitch, uint degSum, GAS_SURF gas)
{
	uint tid = blockIdx.x*blockDim.x + threadIdx.x;
	if(tid>allNum) return;
	uint x = Data[tid];
	uint y = x&0x000007ff;//y=x<<(32-PIX_BIT)>>(32-PIX_BIT);
	x = x>>PIX_BIT;
	if(x>=cRenderParam.renderWidth || y>=cRenderParam.renderHeight)			return;

	float* pixel = (float*)(surface + y*pitch) + 4*x;
	pixel[0] = 1.0;
	pixel[1] = 1.0;
	pixel[2] = 1.0;
	pixel[3] = 0.0f;


	float dir[3] ={2*(x+0.0f)/cRenderParam.renderWidth-1.0f, 1.0f-2*(y+0.0f)/cRenderParam.renderHeight, 0.0f};

	float trans[4];
	trans[0] = dir[0]*cRenderParam.invViewProj[0][0]+dir[1]*cRenderParam.invViewProj[1][0]+dir[2]*cRenderParam.invViewProj[2][0]+cRenderParam.invViewProj[3][0]; 
	trans[1] = dir[0]*cRenderParam.invViewProj[0][1]+dir[1]*cRenderParam.invViewProj[1][1]+dir[2]*cRenderParam.invViewProj[2][1]+cRenderParam.invViewProj[3][1]; 
	trans[2] = dir[0]*cRenderParam.invViewProj[0][2]+dir[1]*cRenderParam.invViewProj[1][2]+dir[2]*cRenderParam.invViewProj[2][2]+cRenderParam.invViewProj[3][2]; 
	trans[3] = dir[0]*cRenderParam.invViewProj[0][3]+dir[1]*cRenderParam.invViewProj[1][3]+dir[2]*cRenderParam.invViewProj[2][3]+cRenderParam.invViewProj[3][3]; 
	dir[0] = trans[0]/trans[3]-cRenderParam.eye[0];
	dir[1] = trans[1]/trans[3]-cRenderParam.eye[1];
	dir[2] = trans[2]/trans[3]-cRenderParam.eye[2];

	float tmin, tmax;
	int firstSeg[3] = {0, 0, 0};
	int seg[3] = {1, 1, 1};
	if(0==RayBoxIntersect(&tmin, &tmax, cRenderParam.eye, dir, gas.first, gas.def, firstSeg, seg)|| tmax<1 ||fabsf(tmin)>1.0e3||fabsf(tmax)>1.0e3)	return;
	if(tmin<1)	tmin=1;

	float fValue[MAX_DEG_INTERPOLATION+1];
#if _KI_BEZ_
	float coeff[MAX_DEG_INTERPOLATION+1];
#else
	float coeff[MAX_CONTROL_POINT];
#endif

	//********* code surround by /**/ is chebyshev interpolation ****************//
	//remember: change the interpolation here will also need change the InvMatrix;
	float In[3] ={cRenderParam.eye[0]+tmin*dir[0],cRenderParam.eye[1]+tmin*dir[1],cRenderParam.eye[2]+tmin*dir[2]};
	float Out[3]={cRenderParam.eye[0]+tmax*dir[0],cRenderParam.eye[1]+tmax*dir[1],cRenderParam.eye[2]+tmax*dir[2]};
	float xyz[3];
	float cheby;
	for(int i=0; i<=degSum; i++)
	{
		cheby = (1-(__cosf(PI*(2*i+1)/(2*(degSum+1)))))/2;
		xyz[0] = (1-cheby)*In[0]+cheby*Out[0];
		xyz[1] = (1-cheby)*In[1]+cheby*Out[1];
		xyz[2] = (1-cheby)*In[2]+cheby*Out[2];
		fValue[i] = EvalGenAlgSurf(xyz[0], xyz[1], xyz[2]);
	}

	uint matrixInd = 0;
	for(int i=0; i<=degSum; i++)
	{
		coeff[i] = 0.0;
		for(int j=0; j<=degSum; j++)
		{
			coeff[i] += fValue[j]*cInvMatrix[matrixInd++];
		}
	}
	int indC;
	for(indC=0; indC<degSum; indC++)
	{
		if(coeff[indC]*coeff[indC+1]<0)	break;
	}

	if(indC==degSum)	return;


#if _KI_BEZ_
 	float firstRoot = KnotInsertBez(coeff, degSum);
#else
 	float firstRoot = KnotInsert(coeff, degSum);
#endif


 	if(firstRoot<0)		return;
	firstRoot = tmin*(1-firstRoot) + tmax*firstRoot;
	xyz[0] = (cRenderParam.eye[0]+firstRoot*dir[0]);
	xyz[1] = (cRenderParam.eye[1]+firstRoot*dir[1]);
	xyz[2] = (cRenderParam.eye[2]+firstRoot*dir[2]);

	float norm[3];
	//float fVal;
	//fVal = EvalGenAlgSurf(xyz[0], xyz[1], xyz[2]);
	//norm[0] = EvalGenAlgSurf(xyz[0]+DERIVATIVE_DELTA, xyz[1], xyz[2]);
	//norm[1] = EvalGenAlgSurf(xyz[0], xyz[1]+DERIVATIVE_DELTA, xyz[2]);
	//norm[2] = EvalGenAlgSurf(xyz[0], xyz[1], xyz[2]+DERIVATIVE_DELTA);
	//norm[0] = (norm[0]-fVal)*DERIVATIVE_DELTA_INV;
	//norm[1] = (norm[1]-fVal)*DERIVATIVE_DELTA_INV;
	//norm[2] = (norm[2]-fVal)*DERIVATIVE_DELTA_INV;
	//norm[0] = -norm[0];
	//norm[1] = -norm[1];
	//norm[2] = -norm[2];
	EvalGenAlgSurfDf(xyz[0], xyz[1], xyz[2], norm);
	normalize(norm);

	pixel[0] = xyz[0];
	pixel[1] = xyz[1];
	pixel[2] = xyz[2];
	pixel[3] = 1.0;
	
	pixel = (float*)(Normal + y*pitch) + 4*x;
	pixel[0] = norm[0];
	pixel[1] = norm[1];
	pixel[2] = norm[2];

	//float rgb[3] = {0.0, 0.0, 0.0};
	////MATERIAL material;
	////chooseMaterial(0, &material, &cRenderParam.material);
	//shading(rgb, xyz, norm, cRenderParam.eye, cRenderParam.nLight, cRenderParam.light, &cRenderParam.material, 0);
	//pixel[0] = rgb[0];
	//pixel[1] = rgb[1];
	//pixel[2] = rgb[2];
}


extern "C" void KIMCOne(RENDER_PARAM* renderParam, unsigned char* surIn, unsigned char* surNorm, unsigned char* surPos, size_t pitch, bool* reNormal)
{
	cudaMemcpyToSymbol(cRenderParam, renderParam, sizeof(RENDER_PARAM));

	dim3 Db = dim3(8,8);
	dim3 Dg = dim3( width/Db.x, height/Db.y);
	if(width %Db.x)	Dg.x++;
	if(height%Db.y) Dg.y++;

	int memS = pixelN*sizeof(uint);	
	CUDA_SAFE_CALL(cudaMemset(dData, 0, memS));
	CUDA_SAFE_CALL(cudaMemset(dValid, 0, memS));	

	//Silhouette<<<Dg, Db>>>(surIn, pitch);
	

	FirstNewton<<<Dg, Db>>>(surIn, surNorm, surPos, pitch, dValid, dData);
	*reNormal = true;


	//uint compNum;
	//compact(&compNum, pixelN*sizeof(uint), dNumber, dData, dValid, dCompact, 0, NULL);  
	////printf("compNum=%u, pixelN=%u\n", compNum, pixelN);

	//uint gridN = compNum/64;
	//if(gridN%64) gridN++;
	//SecondKI<<<gridN,64>>>(compNum, dCompact, surIn, surOut, pitch, degSum, GASData);


}

//extern "C" void KIMCOne(unsigned char* surIn, unsigned char* surOut, size_t pitch)
//{
//	renderParam->material.singleColor = 1;
//	cudaMemcpyToSymbol(cRenderParam, renderParam, sizeof(RENDER_PARAM));
//	
//	//surIn contains the normal info,while surOut position&Sil&contour.
//
//    cudaMalloc((void**)&Debug, 20*sizeof(float));
//
//
//	////Computer 3DDDA
//	dim3 Db = dim3(8,8);
//	dim3 Dg = dim3( width/Db.x, height/Db.y);
//	if(width %Db.x)	Dg.x++;
//	if(height%Db.y)  Dg.y++;
//
//	//dumpBufferfloat((float*)surOut, 512*512*4);
//	//dumpBufferFloatFile((float*)surOut, 512*512*4, 512*4);
//	//dumpBufferfloat(dFillHor, 512*MAX_FILL_N);
//	//dumpBufferfloat(dFillVer, 512*MAX_FILL_N);
//
//	int memS1 = MAX_FILL_N*width*sizeof(float);
//	CUDA_SAFE_CALL(cudaMemset(dFillHor, 0, memS1));
//	int memS2 = MAX_FILL_N*height*sizeof(float);
//	CUDA_SAFE_CALL(cudaMemset(dFillVer, 0, memS2));	
//
//	CUDA_SAFE_CALL(cudaMemset(dNumber, 0, sizeof(size_t)));	
//
//	Classify<<<Dg, Db>>>(surIn, surOut, pitch, Debug, dFillHor, dFillVer, dSilScan, dNumber);
//	//dumpBuffer(dSilScan,hNumber);
//	//dumpBufferfloat(Debug, 16);
//
//	uint hNumber=0;
//	CUDA_SAFE_CALL(cudaMemcpy(&hNumber, dNumber, sizeof(size_t), cudaMemcpyDeviceToHost));
//	//printf("sum=%u\n", hNumber);
//	uint gN = hNumber/64;	if(gN%64)	gN++;
//	UpdataSil<<<gN,64>>>(dSilScan, surIn, surOut, pitch, dFillHor, dFillVer, hNumber, Debug);
//	//dumpBufferfloat(Debug, 20);
//	//dumpBufferfloat(dFillHor, 512*MAX_FILL_N);
//	//dumpBufferfloat(dFillVer, 512*MAX_FILL_N);	
//
//
//	CUDA_SAFE_CALL(cudaMemcpy(hFillHor, dFillHor, memS1, cudaMemcpyDeviceToHost));
//	CUDA_SAFE_CALL(cudaMemcpy(hFillVer, dFillVer, memS2, cudaMemcpyDeviceToHost));
//	CUDA_SAFE_CALL(cudaMemcpy(hImg, surOut, width*height*4*sizeof(float), cudaMemcpyDeviceToHost));
//	float singScreen[2*MAX_SING_PT];
//	GetSingScreen(renderParam, singScreen);
//	SilTopo(singScreen);
//
//
//
//
//	
//	
//	cudaFree(Debug);
//
//}